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Abstract 

We discuss the imaging properties of uniaxial epsilon-near-zero metamaterial slabs with possibly 
tilted optical axis, analyzing their sub-wavelength focusing properties as a function of the design 
parameters. We derive in closed analytical form the associated two-dimensional Greens function 
in terms of special cylindrical functions. For the near- field parameter ranges of interest, we are 
also able to derive a small- argument approximation in terms of simpler analytical functions. Our 
results, validated and calibrated against a full-wave reference solution, expand the analytical tools 
available for computationally-efficient and physically-incisive modeling and design of metamaterial- 
based sub- wavelength imaging systems. 

PACS numbers: 78.20.Ci, 41.20.Jb, 42.30.Wb 
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I. INTRODUCTION 



The seminal work by Pendryi demonstrated that a slab of negative-index material (loss- 
less and impedance-matched with the surrounding medium) would ideally image a perfect 
copy of a source. Such perfect lensing mechanism allows the restoration of the sub- wavelength 
spatial details carried by the evanescent spectral components via efficient coupling to reso- 
nant surface-plasmon states, and may lead to revolutionary applications in a variety of fields, 
including nanolithography, bio-sensing, and spectroscopy At infrared and optical frequen- 
cies, where magnetic activity is difficult to achieve and many materials may naturally exhibit 
negative permittivity, an approximate implementation was suggested 1 and experimentally 
demonstrated^- in the simple form of a thin layer of silver, with severe restrictions on the 
field polarization and near-field range, and the sub-wavelength resolution ultimately limited 
by the material losses. 

The above results have generated a growing interest in the study of metallo-dielectric 



multilayered structures (see, e.g., Refs. l5Ml8l for a sparse sampling), which allow mitigation 



of the above loss- and range-related limitations, as well as further degrees of freedom for 
design optimization. Essentially, these configurations exploit the inherent (possibly extreme) 
anisotropy exhibited by metallo-dielectric multilayers in order to convert evanescent spectral 
components with large transverse wavenumbers into propagating waves. More recently, the 
use of obliquely layered structures has been proposed in order to achieve simple image manip- 
ulation (lateral displacement) with sub- wavelength resolution. 19 Of particular interest is the 
epsilon-near-zero (ENZ) regime, which has also found interesting applications to other sce- 
narios, including cloaking,— >2i light funneling through sub- wavelength apertures,— controlled 
leaky-wave radiation, 23 suppression of Anderson localization in disordered multilayers,— loss- 
induced omnidirectional bending,— and nonlinearity enhancement.—^ Also worth of men- 
tion are the studies on sub- wavelength imaging devices based on wire media (see , e.g., Refs. 



28 



31 



33|). 



30), as well as on the transformation-optics paradigm (see, e.g., Refs. 

With the exception of few cases (see, e.g., Refs. Ls| andll7|), for which analytical approxi- 
mations of the Green's function can be worked out, the imaging properties (e.g., resolution) 
of the above configurations need to be assessed numerically. In this paper, we study in 
more detail this geometry and we show that the two-dimensional (2-D) Green's function 
of a slab of uniaxial ENZ metamaterial with tilted optical axis can be calculated analyti- 



cally in closed-form, greatly facilitating the analysis and design of its anomalous imaging 
properties. Our results are expressed in terms of special cylindrical functions that can be 
efficiently computed, and are also amenable to simple approximations in the parameter 
range of interest. These findings shed new light on the sub-wavelength imaging properties of 
ENZ metamaterial slabs and allow tailoring their properties without the need of extensive 
numerical simulations. 

The rest of the paper is organized as follows. In Sec. [Til we introduce the problem 
geometry and formulation. In Sec. IHH we derive the analytical solution for the Green's 
function (with the more involved details relegated in Appendices IAT 4C]) , discuss the related 
computational issues, and work out more manageable approximations for specific values of 
the design parameters of interest. In Sec. IIVI we illustrate some representative results in 
order to validate and calibrate our proposed solutions, and discuss relevant physical insights 
on the imaging properties of the considered metamaterial lenses. Finally, in Sec. |Vl we 
provide some brief conclusions. 



II. PROBLEM GEOMETRY AND FORMULATION 



A. Generalities 



Figure [T] shows the geometry of the problem. We consider a 2-D scenario featuring a 
metamaterial slab of thickness d, infinitely extent in the y, z directions, immersed in vacuum. 
The metamaterial is assumed to be nonmagnetic and uniaxially-anisotropic with optical axis 
tilted of an angle a with respect to the x-axis. Accordingly, the corresponding permittivity 
tensor is most naturally described in the rotated (principal) reference system (£, v, z) 



e e 
e v 
e v 



(1) 



and represents a rather general homogenized (effective-medium-theory) model for diverse 
classes of artificial materials, such as multilayered structures or wire media. In what fol- 
lows, we assume that the homogenized constitutive parameters in ([T]) do not depend on the 
wavevector. While acknowledging the implied limitations in predicting nonlocal effects that 
can take place in metallo-dielectric multilayers (see, e.g., Refs. I34M36I) or wire media (see, 



e.g., Refs. I37H39I). we focus here on this simplified model which is amenable to analytical 
treatment. 

The essential kinematical (wavevector, group velocity) properties of wave propagation in 
such a medium may be qualitatively understood by looking at the equi-frequency contours 
(EFCs). Once again, these are most easily studied in the spectral- wavenumber plane (kg, k v ) 
associated with the rotated reference coordinate systems (£, v) in Fig. (TJ where they assume 
the canonical form 

k 2 b 2 

f + f = kl (2) 

with k = Uy/Eofio = 27r/Ao denoting the vacuum wavenumber (and Ao the corresponding 
wavelength). At variance with the circular shape exhibited by isotropic media, depending on 
the sign of and e v , these EFCs may be either elliptic (e^e v > 0) or hyperbolic (e^e v < 0), 
as exemplified in Fig. [2j This implies, especially in the hyperbolic case [cf. Fig. EJ^b)], 
that spectral components characterized by large transverse wavevenumbers (which would 
be otherwise evanescent) may actually propagate. Moreover, in the ENZ limit (e v — > 0) of 
interest, the EFCs tend to become much flatter along the k v direction, thereby implying that 
the allowed group velocities (normal to the EFCs) tend to be directed along the £ direction, 
so that sub-wavelength details can be transported along the optical-axis direction. Such 
condition (in conjunction with — > oo), which was first proposed in Ref. [5| in order to 
mitigate the loss-induced limitations in single-layer silver superlenses, may be attained for 
multilayered structures from standard mixing formulas^ using constituent materials with 
opposite-signed permittivities (e.g., noble metals and dielectrics, at optical frequencies). 



This condition was also exploited in Ref. 
The reader is referred to Refs. 



10 



11 



and 



for far-field sub-wavelength imaging. 
12] for different extreme-anisotropy-based sub- 



wavelength imaging systems relying on Fabry-Perot resonance effects. 



B. Green's Function 

From the above observations, we expect that the uniaxial ENZ metamaterial slab of in- 
terest may exhibit sub-wavelength image formation and lateral-displacement capabilities. 
In what follows, we analytically study the response to a unit-amplitude, time-harmonic 
[exp(— iuit)), z— directed magnetic-current (V/m 2 ) line source located at a distance x s (typ- 
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ically <C Aq) from the slab (cf. Fig. [TJ) 



M z (x,y) = S (x + x s )S(y) 



(3) 



We start recalling the well-known Green's function of vacuum for the z— directed magnetic 
field, 

Gf\x,y;x s ) ■■ 



w£o H (i) 
4 



koV( x + x s ) 2 + y 2 



ujeq f exp [i {k x \x + x s \ + k y y)} 



An 



dk 



(4a) 
(4b) 



in terms of the zeroth-order Hankel function of first kind (cf. Sec. 9 in 41), or the 
corresponding spectral-integral representation, 42 with 



k x = Jk 2 -k 2 , lm(^)>0. 



(5) 



Accordingly, the presence of the metamaterial slab can be accounted for by introducing in 
the spectral integral representation (14bp the corresponding transverse-magnetic plane-wave 
transmission-coefficient T(k y ) and an appropriate displacement along the x direction, viz., 



4n 



exp [i (k x \x + x s — d\ + k y y)] dk y . 



(6) 



The spectral integral in cannot be generally calculated analytically in closed form. In 
certain regimes (see, e.g., Refs. |8| and Il7l). closed- form near- field approximations may be 
worked out by applying Cauchy's residue theorem and neglecting the branch-cut contribu- 
tion. However, in general, a brute-force numerical integration is needed. 



C. ENZ Regime 

In the limit e v — > of interest, the transmission coefficient to be considered in ([6]) reduces 
to (see Appendix |A] for details) 

T n \__ 2ggfcxexp (-ikydtana) 
y 2e^k x — ide^k 2 + idk 2 sec 2 a' 

from which it can be observed that total-transmission (i.e., perfect impedance matching) 
can only be achieved for 

k y = ±kQ^/el cos a. (8) 
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The above condition can be satisfied for propagating waves impinging from vacuum (i.e., 
\k y \ < ho) with incidence angle (see Fig. [[]) 



9 = ± arcsin (,,/e£Cos a) , 

which admits real solutions for positive values of and .y/££ |cosa| < 1. 

The transmission coefficient in ([7]) can be recast in a convenient canonical form 

Ak x exp (—ikyd tan a) 



(9) 



T(ky) 



where 



A 



[k x - ki) {k x - «a) 

2ie^ cos 2 a 
= d ' 



«1,2 = - 



cos 2 a 



d 



(10) 
(11a) 

;iib) 



ie^ =F y k\d 2 sec 2 a (sec 2 a — e^) — e 2 - 

Accordingly, by substituting ( ITOj) into (jSJ), the slab response is reduced to calculating the 
canonical spectral integral 

oo 

Aue f exp [i (k x x d + k y y d )} 



Air 



(k x - Kx) (k x - « 2 ) 



■dk 



with (15]) and 



= x + x s — d, i)d = y — d tan a. 



(12) 



(13) 



III. ANALYTICAL RESULTS 



General Solution 



Introducing the (normalized) polar coordinates 



( = k Jx 2 d + y% tp = arctan — 



Ud 



(14) 



it can be shown (see Appendix [B] for details) that the canonical spectral integral in f fT2|) 
admits closed-form solutions of the type 

Acu^o 



Gf\x,y- 1 x s ) 



47r(/«i - K 2 ] 



(15) 



where 



FiAC, f) = exp (aJaC) ^2 + Xi,2^-% ( a i^ Vi,2, C) 
+ exp (a^ 2 C) ^r,2 + Xi,2H#o { a i,2, Vi,2, C) 



(16) 
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In ( TT6|) . A± 2 and A 12 are integration constants, 



i ( K i,2 cos ip ± I sin if \ \Jk$ — k\, 



a% = ^- ^, Im ( ^/*g - < 2 ) > 0, ( 1 7) 



71 Ki 2 ( COS f\ kn — k\ 2 =F | Sill 99 1 Ki 2 

*£. = — ^ V ' (18) 



are known parameters depending on the frequency, slab parameters, and observation direc- 



tion, and 



Hek 1) (a,ri,0= f exp (-ar) ^\r)dr (19) 
Jn 

denotes the complementary incomplete Lipschitz-Hankel integral (CILHI) of the Hankel 
type^ - a special function belonging to the general class of incomplete cylindrical function.— 
These functions have insofar found applications in other areas of electromagnetics, ranging 
from diffraction to traveling- wave sources (see, e.g., Ref. |45| for a review) 

The lower limit of integration 77 in ffl9l) is chosen so as to guarantee convergence at infinity 
in the complex r plane, viz., 

J 00, Re(a) > 0, 

I ooexp(z7r), Re(a) < 0. 

The calculation of the 92-dependent integration constant A± 2 and A{ 2 in (fl6l) is generally 
quite involved (see the discussion in Appendix [C]) . Nevertheless, for the near-field config- 
uration of direct interest to this investigation, featuring a source located very close to the 
input slab interface (x s <C Ao), this calculation becomes remarkably simpler for two specific 
observation planes located in the two principal planes of the image, namely: transverse cuts 
at the output slab interface (x = d), and orthogonal cuts passing through the fiducial image 
position (y a = rftana). For the former case (x — d), we obtain (see AppendixIClfor details) 

A{ 2 = 0, (21a) 

A i,2 = Ci,2 - Xx^o i a t,2> Vij, °) - Xi^^o ( fl 1^2> Vl,2, °) > ( 21b ) 



where 



L\ 2 = ; : < 2arctan 



lh.2 _ K 2 \ U2 _ 2 

tv h, 12 ^ \y o 1,2 / 

+ 7r{l-4«[Re(« li2 )]M[Im(K lj2 )]} I, (22) 



with u denoting the Heaviside unit-step function, ano~ 



v (!) ( n\ 2ilog ( Q + V q2 + !) ~* (0 o\ 
Heb'(a,T},0) = -==== , (23) 

7TVGr + 1 



with the principal branch chosen for the natural logarithm, and the branch-cut for -y/a 2 + 1 
chosen as 

Re (y/aFTT) > 0, Va, 

~ (24) 



Im (Va 2 + 1) > 0, for Re (y/cF+T) = 0. 

For the latter case (y = dtana), the above expressions are still valid if Re(a^ 2 ) = Re(a]~ 2 ) < 
0, otherwise they simply reduce to (see Appendix O for details) 

A+ 2 = A~ lt2 = 0. (25) 

It is rather remarkable that in these two relevant planes, of most interest to analyze the 
imaging properties of the metamaterial slabs, we are able to obtain a closed-form analytical 
solution of the field distribution induced by an arbitrary magnetic source. This result will be 
used in the following to highlight the remarkable imaging properties of ENZ metamaterials 
with arbitrarily tilted anisotropy axis. 



B. Physical Interpretation 

It is of particular interest to identify the distinct phenomena involved in the wave inter- 
action with the geometry of Fig. [1] as a function of the various parameters considered in our 
analytical solution in ( TlBl) and ( Tl6|) . First, we notice that the exponential terms exp (af 2 C) 
account for the dominant resonances exhibited by the slab, corresponding to pole singu- 
larities in the spectral integral formulation.— Their localization properties can be readily 
related to the slab physical parameters via ffTTl) and (lllbl) . 

The terms exp (af 2 () (aj 2 , r}f 2 i C) instead account for the continuous radiation 
spectrum, corresponding to the branch-cut contributions in the spectral integral formulation.— 
We remark that, unlike the configuration in Refs. [si and 17, these contributions are generally 
non- negligible. Recalling the definition in ( fl9l) . these terms may be physically interpreted as 
smoothed versions (via convolution with a complex-exponential window) of a zeroth-order 



Hankel function of the first kind H 



(i) 



C. Computational Aspects 



1. Numerical Computation of CILHIs 

Recalling the expression of the af 2 parameters in (fT7|) . it is evident that the numerical 
implementation of our analytical solution requires in general the calculation of complex- 



argument CILHIs of the Hankel type. In Ref. 



43 



various series expansions were derived for 



accurate and efficient computation of these special functions. In particular, by comparison 
against brute-force numerical integration (via adaptive Gauss/Kronrod quadrature) of the 
corresponding spectral integrals, efficiency improvements ranging from one to nearly three 
orders of magnitudes were found. 

Our numerical implementation is based on a selective application of the factorial 
Neumann and Struve-function series expansions, following the guidelines of Ref. 



43 



for 



the various parameter ranges. Although a proper numerical optimization to maximize the 
calculation efficiency of these coefficients goes beyond our interest, it is evident that our 
closed-form analytical solution may significantly outperform conventional numerical solvers. 



2. Small- Argument Approximation 



Since the application of the metamaterial slab and the focus of this paper is concentrated 
on the near-field sub-wavelength imaging scenario, with x s <C Ao and x = d, a simple small- 
argument (£ <C 1) analytical approximation for the CILHIs may be conveniently utilized. 
First, we recall the small- argument approximation of the zeroth-order Hankel function of 
the first kind (cf. Eqs. (9.1.12) and (9.1.13) in Ref. M), 



H? } (r) ~ 1 + ( ~- 



2i 



71 



!og ( g ) + 7 



(26) 



where 7 denotes the Euler-Mascheroni constant.— Using this approximation, it follows that 



exp (aQ I exp (— ar) Hq 1 ^ (r) dr 
'0 



exp (aQ — 1 



Tia 6 V2 



2i {exp «) [E, «) + log (2a)] + 7} 



7ra 



(27) 



where E\ denotes the exponential integral [cf. Eq. (5.1.1) in Ref. |4l| which, recalling the 



expansion in Eq. (5.1.11) of Ref. |4l| and exploiting a (1,1) Pade approximant^ for the 



power-series part, may be conveniently approximated as 



E 1 «) 



4< 



- 7 - log «) 



4 + aC 

Recalling ( II 9p and ( I23p . we finally obtain the small-argument approximation 

2i log (a + Va 2 + l) - vr 8iaC 



exp(aC)'He[ ) 1) (a, rj, () ~ exp(aC) 



[exp «) — 1] < 2z 



vrVa 2 + 1 

c 



7ra (4 + aC) 



log 



+ 7 



+ 7T 



7ra 



in terms of simple analytical functions. By substituting ( 129|) in ( IT6|) . we obtain 



F lj2 ((,<p) ~ A+ 2 exp(a+ 2 C) +xi 

"2*log(a+ 2 + x /(a+ 2 ) 2 + l 

+ x7 2 ex P( a i!2C) 



[exp (a^ 2 C) - l] < 


[• 


log 


(9- 


+ 7r| 


7ra] 


\- 

.2 



7T 



8ia^ 2 C 



W« 2 ) 2 + l 



7ra+ 2 (4 + a+ 2 C) 



[exp (a 12 C) - 1] S 2i 



+ X 



log 



+ 7 



7T 



1,2 



ixa 



1,2 



+ Xl,2 eX P( a l,2C) 



2ilog f a 12 + W(a lj2 ) 2 + 1 



7T 



8ia 12 C 



1,2 ( 4 + a l, 2 C) 



(28) 



(29) 



(30) 



which, substituted in (fT5l) . yield the final approximation for the Green's function (not given 
explicitly here for brevity). Once implemented, this solution may provide a complete de- 
scription of the imaging properties of the metamaterial slab, based on conventional basic 
functions. Its overall applicability and accuracy, which is expected to be restricted to small 
values of the argument, will be quantitatively assessed in Sec. IIV Bl 



IV. REPRESENTATIVE NUMERICAL RESULTS AND PHYSICAL INSIGHTS 
A. Generalities 

In what follows, for certain representative ENZ parameter configurations, we validate and 
calibrate the analytical solutions derived in Sec. II I II against a reference solution obtained via 
brute-force numerical quadrature of the spectral integral in ([6]) with the general transmission 
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coefficient in ( lAlj) . We further discuss the imaging properties of the metamaterial slab 
for specific design parameters of interest. Our reference solution relies on Gaussian- type 
quadrature rules (cf. Sec. 25.4 in Ref. 4lL with the number of nodes refined adaptively so 
as to guarantee a four-digit accuracy. It is worth pointing out that no particular attempt was 
made to optimize its numerical implementation, since a careful and thorough assessment of 
the computational convenience of CILHI-based solutions vs. numerical spectral integration 
was already carried out in Ref. 43| (see also the discussion in Sec. IIII C lj) . and is therefore 
not the focus of our investigation here. 



B. Numerical Results and Discussion 

We start considering an ideal lossless configuration featuring a line source placed at a 
distance x s = Ao/50 away from a slab of thickness d = 0.5Ao, and constitutive parameters 
e v = 0, £g = —2. Figure [3] shows the normalized Green's function intensity maps, within 
and beyond the image plane x = d, computed by using the reference solution [cf. (|5J], 
for values of the optical-axis angle a ranging from to 75°. Similarly to what observed in 
Ref. |l9|, the images are maximally localized at the image plane x = d and are subject to a 
lateral displacement (of a quantity y a = dt&na) with respect to the source position (y = 0), 
which is in agreement with the fact that sub-wavelength details of the source are effectively 
transported with small distortion along the optical-axis direction. As it can be observed, 
increasing values of the tilt angle a (i.e., of the lateral displacement) are accompanied by 
progressively worse localization properties (note the different color scales of the plots). 

Figures H] and [5] show the intensity profiles along the two principal cuts, namely, the 
image plane at the slab interface (x = d) and the orthogonal plane passing trough the 
fiducial image position (y a = rftana), respectively, comparing the reference solution [cf. 
([6]), black-solid curves] with the CILHI-based analytical solutions in (jT5l) and (|T6|) [with the 
parameters A and given by pip ]. Along these cuts, it was possible to apply the small- 
argument approximation and derive Eqs. (1291) . More specifically, both the "exact" (i.e., 
numerically computed as in Sec. IIII C II) and small-argument-approximated [cf. ( 1291) ] CILHIs 
are considered (red-dashed and blue-dotted curves, respectively). The "exact" CILHI-based 
solutions are practically indistinguishable from the reference solutions [cf. ([6])] on the scale 
of the plots, whereas, as expected, the small-argument approximation breaks down away 
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from the peak. These limitations are more evident in the orthogonal cuts (Fig. [5]), since the 
longitudinal localization is typically poorer. It is important to stress that, compared to a 
conventional focusing lens, the metamaterial slab is able to focus the transverse image well 
below the diffraction limit, but it is not very effective in focusing in the longitudinal plane. 
Effectively, as shown in Fig. [5J the sub-wavelength spot decays away from the slab with 
a conventional exponential drop, due to the diffraction of the high-wavenumber spectral 
components of the image transported to the back of the slab. 

It is interesting to notice in Figs. H] and how the small- argument approximation is 
sufficient to correctly capture the essential localization properties of the image, especially 
in the transverse direction. In what follows, for quantitative assessments, we consider two 
typical figures of merit: the f ull- width- at-half -maximum (FWHM) and the (normalized) 
peak intensity at the image plane x = d. Figure |6] compares the FWHM and (normalized) 
peak- intensity, estimated via the reference solution [cf. (|6]), circular markers] and the small- 
argument CILHI-based analytical solution (square markers), as a function of a. Overall, 
a uniformly good agreement is observed, with maximum errors < 2% in the FWHM and 
< 0.3% in the peak-intensity. Consistently with the visual impression from Figs. |3]andHl 
both observables deteriorate with increasing values of a. More specifically, over the range 
< a < 75°, the FWHM increases from ~ 0.12Ao to ~ 0.35Ao, while the (normalized) peak- 
intensity decreases from ~ 3.3 to ~ 0.2. To sum up, the ENZ configuration analyzed in this 
paper is able to transport, and laterally displace, sub-wavelength details of a source from 
the input to the output slab interface, with resolutions as good as a tenth of a wavelength. 
However, large lateral displacements result in a larger degradation of the image resolution 
and intensity. For instance, lateral displacements of nearly two wavelengths [cf. Figs. [3]T), 
H]T), and[5]T)] may be attained at the expense of a factor ~ 3 in the resolution and over an 
order of magnitude in the intensity. 

For the same configuration, and a fixed optical-axis direction (a = 0), Fig. [7] shows the 
two considered figures of merit as a function of in the hyperbolic-medium regime. As 
expected, recalling our discussion at the end of Sec. Ill Al the figures of merit improve for 
larger absolute values of e%. Conversely, they strongly deteriorate for e% — > 0. This is also 
not surprising, as it is well known that in the isotropic ENZ limit the slab becomes a highly- 
selective spatial filter.— Also in these cases, the agreement between the reference solution 
[cf. ()6])] and the small- argument CILHI-based analytical solution is very good, with same 
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maximum errors as above. 

Next, we assess the critical ENZ assumption underlying our analytical solution. It is well 
known that in practical implementations (e.g., metallo-dielectric multilayers) the condition 
e v = may only be approximatively attained due to the presence of losses, and in any 
case limited to one single frequency point. While it is possible, in principle, to achieve 
a vanishing real part (at a given frequency), zeroing the imaginary part is prevented by 
unavoidable material losses. Nevertheless, recent studies^ have demonstrated the promising 
potentials of gain-assisted nanocomposites for the synthesis of artificial materials with very 
small values (~ 10~ 2 ) of the real and imaginary part of e v . In order to assess the applicability 
of our proposed analytical solution to such regime, we consider a more realistic configuration 
with, d = 0.5Ao, ct = 0, Re(e^) = —2, and a small but nonzero \e v \. More specifically, we 
assume Re(e v ) = 10~ 3 , and let the imaginary parts of e% and e v vary over several decades. 
Figure [8] shows the corresponding FWHM and peak-intensity estimated via the reference 
solution [cf. ()6])]. As expected, for asymptotically vanishing losses, they approach the 
estimates from our small-argument CILHI-based analytical solution (shown as horizontal 
dashed lines), and progressively depart from them for increased loss levels. In particular, for 
lm(e^ v ) = 10 -2 , the agreement is still satisfactory, with only a ~ 2% error in the FWHM 
and a ~ 3% error in the peak-intensity. For lm(e^ v ) = 0.05, the errors increase to ~ 15% 
and ~ 16%, respectively, which may still be acceptable. Qualitatively similar trends where 
observed for different values of Re(e^), as exemplified in Fig. [9l In this case, pertaining 
to Re(e^) = —5, the errors in the FWHM and peak-intensity are ~ 0.4% and ~ 1%, 
respectively, for Im(e^u) = 10~ 2 , and ~ 5% and ~ 9%, respectively, for lm(e^ v ) = 0.05. 
From the physical point of view, our results indicate that the imaging properties of the ENZ 
anisotropic slabs considered here are quite robust to losses and frequency variations. 

Overall, the above results also indicate that our proposed analytical solution accurately 
captures the essential image- formation properties in the scenario of interest, and it can be 
safely applied in low-loss scenarios. 

V. CONCLUSIONS 

In this paper, we have presented an analytical study of the sub-wavelength imaging 
properties of uniaxially-anisotropic ENZ metamaterial slabs with tilted optical axis. In 
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particular, we have derived a closed-form analytical solution for the 2-D Green's function 
in terms of special cylindrical functions. These functions can be efficiently computed via 
well-established numerical schemes, yielding computational savings up to nearly three orders 
of magnitudes by comparison with brute-force numerical quadrature of the corresponding 
spectral integrals. Moreover, in the near-field parameter range of interest, they can be 
conveniently approximated in terms of simple analytical functions. 

Validation and calibration of our results against a numerical-integration-based reference 
solution [cf. (jBJ)] confirmed the applicability of our solution to sub-wavelength imaging 
scenarios with low-loss constitutive parameters that are within reach of current (e.g., gain- 
assisted) technologies. We have employed this solution to analyze the imaging properties of 
anisotropic ENZ metamaterial slabs varying the design parameters and tilt angle. 

Current and future investigations are aimed at exploiting our proposed parameterization 
in design procedures and optimization schemes to improve the imaging capabilities of these 
devices. Also of interest are possible extensions to account for spatial-dispersion (nonlocal) 
effects, as well as the to predict the emission enhancement for quantum sources radiating in 



Appendix A: Pertaining to (FH) 

The plane- wave transmission coefficient of the uniaxial slab described by the permittivity 
tensor in ([T]) is computed by expanding the field inside the slab in terms of forward and back- 
ward plane waves (with conserved transverse wavevectors) , and enforcing the phase- matching 
and transverse-field-continuity conditions at the interfaces. For the assumed transverse- 
magnetic polarization, via cumbersome but straightforward algebra, we obtain the general 
expression 
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represent the Cartesian components of the relative permittivity tensor in (pQ), and 



£xyk y =F \ i^xx^yy £ X y) { £ xx^q /c 2 ) 

k x i,2 = , (A3) 

i^y) = ^xy^x &xx i^yykx i k X(; ) ^ £ X yky, £ = 1, 2. (A4) 

It can readily be verified that the limit e v — > leads to a 0/0 indeterminate form in (IA1|) . 
In order to evaluate this limit, noting from fl A3 j) that k x \ — > k x2 , it is expedient to linearize 
the exponential function in the denominator of (lAlj) . viz., 

1 yj W C/r(A: y )f/ 2 + (A: y ) - [1 + z - M rf] U+(k y )U 2 (k y ) ' 1 j 

Next, via straightforward McLaurin expansions, we obtain for the various terms in (IA5I) 



£ ra - £xx£w) (kxi - k x2 ) ~ 2 (e^) » ^/ fcgej cos 2 a-k 2 + O U% V (A6) 



Ur(ky)U+(k y ) - U+{k y )U 2 {k y ) ~ 4 foe*)* A^A; 2 ^ cos 2 a - A; 2 + O (4) , (A7) 
(A;^ - fe x i) U+(k y )U 2 (k y ) ~ 2^4 (*ge £ cos 2 a - A; 2 ) 1 sec 2 a + O (e 2 ) , (A8) 
with O(-) denoting the Landau symbol. The limit in (JTj) readily follows by substituting 

3 

(1A6[) f]A8h into fjA5[) . neglecting higher-order terms, and simplifying the dominant term el- 



Appendix B: Derivation of the General Solution in (115ft — (I18j) 



First, by rewriting the rational part of the integrand in (JT2"|) as 

1 «4 K 2 



(k x - Ki)(k x - k 2 ) k x (k x - /ci)(«i - k 2 ) k x (k x - k 2 )(k 1 - k 2 ) 
we obtain [cf. f fl5|) ] a different, generic canonical integral of the form 



exp 



iQ — cos tp + -f- sm <z> 

fc ft 



k x {k x k) 
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(Bl) 



(B2) 



5CH . closed-form calculation of spectral 



in the polar coordinates given by (JHJ). In Ref. 
integrals of this type was carried out via rather cumbersome contour-integration techniques. 
In what follows, we rely on an alternative, relatively simpler, differential-equation-based 



procedure proposed in Ref. 



51 



(based on the work in Ref. |52|). The basic idea is to 
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construct an inhomogeneous differential equation (in the ( variable) satisfied by the integral 
in (IB2j) . To this aim, we apply to (IB2j) a second-order differential operator 



V 2 [F] (C, <p) = 



a 



a 



ac 2+/?1 9C +/3 ° 

K cos(2<^) 



(B3a) 



k 2 



cos tp + k y sin y?) 

K 



o /c x A;„ sin(2o9) 
sin 



/.• 2 



(B3b) 



where /?o and 0% are ^-dependent coefficients, and (lB3bl) follows from straightforward differ- 
entiation under the integral sign in flB2j) . Particularly expedient is to choose the coefficients 

as 



fa = sin 2 (p- -3, /9i 
fc 



2w cos 



which allows recasting ( lB3h in a simplified form 



^ rwi/*. \ f« — cos(2o?) — fc„ sin(2<z>)l (k x — k) _ 
V 2 [F] (C, <p) = [ - 1 ^ * 1 rnK x J -F(£, <p) 



k 2 



(B4) 



(B5a) 



k\k — k x cos(2y9) — k y sin(2<^)] 



x 



cxp 



k 2 

«C — cos ip + — - sin ip 
k k 



(B5b) 



Recalling the spectral-integral representation of the zeroth-order Hankel function of the 
first kind [cf. (|4b|)]. 



exp 



«C i— cos if + — ^ sin ip 
k k 



dk 



(B6) 



the right hand side of (IB5I) can be readily calculated, yielding 











(B7) 



We have thus demonstrated that the canonical integral in flB2p may be alternatively cal- 
culated by solving the second-order, inhomogeneous differential equation in (IB7P [with D 2 
given in (1B3al) with (1B4j) ]. This can be accomplished in a systematic fashion by applying 
the method of variation of parameters,™ yielding 



F(C, ip) = A + exp (a + C) + A~ exp (a"C) + Q + ((, <p) + Q"(C, ¥>), 



(B8) 
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where A + and A are integration constants to be determined (see Appendix ICl below) . 



fff 



are the roots of the characteristic equation, and 



Q ((, <p) 



1 p 

— y ± exp [a ± (C - r)] F (r, <p)dr 



iTTKCOSV? .Jj(l)/^ 



+ 



7r/t (k + ikoa cos 



exp (a ± C) fte^ (aW,C) 



(B9) 



(BlOa) 



(BlOb) 



represents the particular solution, with the second equality following from (1B7|) [recalling 

(TO- 

The general solution in (I15p — (j!8p immediately follows, by recalling (IBip and rearranging 
terms.— 



Appendix C: Calculation of the Integration Constants in (1211) 

The ^-dependent integration constants Af 2 in fT2T|) can be determined by enforcing the 
proper boundary conditions (usually, at ( = and for £ — > oo). The general calculation 
procedure is rather involved, and depends on the sign of the real parts of the af 2 parameters 
in ( ITT]) . For instance, if these real parts are all negative, the conditions for ( — > oo cannot 
be exploited, and one is led to enforce only the boundary conditions at £ = 0. However, 
the calculation becomes particularly cumbersome, as the functions F 12 in f)16p are nondif- 
ferentiable at £ = 0. Nevertheless, for the near-field scenario of interest here, the calculation 
becomes particularly simple in the two principal planes of the image, namely x = d and 
y = dtana, which are the most interesting in order to assess the transverse and longitudinal 
localization. 

More specifically, assuming x s = and x = d, we note that, from flTTJ), 

Re (a^ 2 ) = -Re (a+ 2 ) > 0, (CI) 
Accordingly, (12 lap readily follows by enforcing the regularity-at-infinity conditions 

lim F li2 (C,^) = 0. (C2) 
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The remaining integration constants are computed by enforcing the conditions at ( = 0, 
viz., 

F li2 (0, if) = A+ 2 + xt^e^ « 2 , Tjt t2 , 0) + Xi,2^4 ] fe, %2, 0) (G3a) 

oo 

= «i,2 / 7-77 7^2/ = Ci,2, (C3b) 



where the second equality follows from ( 1B2j) . and the arising spectral integral admits the 
closed- form analytical solution given in fl22|) . The above derivations also hold approximately 
for source positions very close to the input slab interface (x s <C Ao). 

For observations along orthogonal planes passing through the fiducial image position 
(yd = 0, i.e., y = tana), we note from (TT7T) that 

a i,2 = ( C4 ) 

which basically means that there are only two effective integration constants to be deter- 
mined. Accordingly, if Re (af 2 ) < 0, we can arbitrarily set A{ 2 = and proceed as above 
[cf. ( 1C3I) ]. so that the results in ( 12TI) hold for this case too. Otherwise, if Re (af 2 ) > 0, then 
the results in (1251) readily follow from the regularity- at- infinity conditions. 



* [ygaldi@uni sannio.it 

1 J. B. Pendry, Phys. Rev. Lett. 85, 3966 (2000). 

2 Z. Liu, N. Fang, T.-J. Yen, and X. Zhang, Appl. Phys. Lett. 83, 5184 (2003). 

3 D. Melville and R. Blaikie, Opt. Express 13, 2127 (2005). 

4 N. Fang, H. Lee, C. Sun, and X. Zhang, Science 308, 534 (2005). 

5 S. A. Ramakrishna, J. B. Pendry, M. C. K. Wiltshire, and W. J. Stewart, J. Mod. Opt. 50, 
1419 (2003). 

6 S. A. Ramakrishna and J. B. Pendry, Phys. Rev. B 67, 201101(R) (2003). 

7 W. Cai, D. A. Genov, and V. M. Shalaev, Phys. Rev. B 72, 193101 (2005). 

8 B. Wood, J. B. Pendry, and D. P. Tsai, Phys. Rev. B 74, 115116 (2006). 

9 K. J. Webb and M. Yang, Opt. Lett. 31, 2130 (2006). 

A. Salandrino and N. Engheta, Phys. Rev. B 74, 075103 (2006). 

1 P. A. Belov and Y. Hao, Phys. Rev. B 73, 113110 (2006). 

18 



12 X. Li, S. He, and Y. Jin, Phys. Rev. B 75, 045103 (2007). 

13 P. A. Belov, Y. Zhao, Y. Hao, and C. Parini, Opt. Lett. 34, 527 (2009). 

14 B. Zeng, X. Yang, C. Wang, Q. Feng, and X. Luo, J. Opt. 12, 035104 (2010). 

15 R. Kotynski and T. Stefaniuk, Opt. Lett. 35, 1133 (2010). 

16 Y. Jin, Prog. Electromagn. Res. (PIER) 105, 347 (2010). 

17 H. Liu and K. J. Webb, Opt. Lett. 35, 1869 (2010). 

18 R. Kotynski, T. Stefaniuk, and A. Pastuszczak, Appl. Phys. A-Mater. 103, 905 (2011). 

19 J. Wang, H. Y. Dong, K. H. Fung, T. J. Cui, and N. X. Fang, Opt. Express 19, 16809 (2011). 

20 A. Alu and N. Engheta, Phys. Rev. E 72, 016623 (2005). 

21 D. S. Filonov, A. P. Slobozhanyuk, P. A. Belov, and Y. S. Kivshar, Phys. Status Solidi 6, 46 
(2012). 

22 D. C. Adams, S. Inampudi, T. Ribaudo, D. Slocum, S. Vangala, N. A. Kuhta, W. D. Goodhue, 
V. A. Podolskiy, and D. Wasserman, Phys. Rev. Lett. 107, 133901 (2011). 

23 K. Halterman, S. Feng, and V. C. Nguyen, Phys. Rev. B 84, 075162 (2011). 

24 A. A. Asatryan, L. C. Botten, M. A. Byrne, V. D. Freilikher, S. A. Gredeskul, I. V. Shadrivov, 
R. C. McPhedran, and Y. S. Kivshar, Phys. Rev. B 85, 045122 (2012). 

25 S. Feng, Phys. Rev. Lett. 108, 193904 (2012). 

26 C. Argyropoulos, P. Y. Chen, G. DAguanno, N. Engheta, and A. Alu, Phys. Rev. B 85, 045129 
(2012). 

27 C. Argyropoulos, P. Y. Chen, F. Monticone, G. DAguanno, and A. Alii, Phys. Rev. Lett. 108, 
263905 (2012). 

28 P. A. Belov, C. R. Simovski, and P. Ikonen, Phys. Rev. B 71, 193105 (2005). 

29 P. A. Belov, Y. Hao, and S. Sudhakaran, Phys. Rev. B 73, 033108 (2006). 

30 P. A. Belov and M. G. Silveirinha, Phys. Rev. E 73, 056607 (2006). 

31 M. Tsang and D. Psaltis, Phys. Rev. B 77, 035122 (2008). 

32 W. Yan, M. Yan, and M. Qiu, J. Opt. Soc. Am. B 26, B39 (2009). 

33 I. Gallina, G. Castaldi, V. Galdi, A. Alu, and N. Engheta, Phys. Rev. B 81, 125124 (2010). 

34 J. Elser, V. A. Podolskiy, I. Salakhutdinov, and I. Avrutsky, Appl. Phys. Lett. 90, 191109 
(2007). 

35 A. A. Orlov, P. M. Voroshilov, P. A. Belov, and Y. S. Kivshar, Phys. Rev. B 84, 045424 (2011). 

36 G. Castaldi, V. Galdi, A. Alu, and N. Engheta, Phys. Rev. Lett. 108, 063902 (2012). 

19 



37 P. A. Belov, R. Marques, S. I. Maslovski, I. S. Nefedov, M. Silveirinha, C. R. Simovski, and S. 
A. Tretyakov, Phys. Rev. B 67, 113103 (2003). 

38 C. R. Simovski and P. A. Belov, Phys. Rev. E 70, 046616 (2004). 

39 R. J. Pollard, A. Murphy, W. R. Hendren, P. R. Evans, R. Atkinson, G. A. Wurtz, and A. V. 
Zayats, V. A. Podolskiy, Phys. Rev. Lett. 102, 127405 (2009). 

40 A. Sihvola, Electromagnetic Mixing Formulas and Applications (IEE Publishing, London, UK, 
1999). 

41 M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Ninth printing (Dover, 
New York, 1970). 

42 L. B. Felsen and N. Marcuvitz, Radiation and Scattering of Waves (IEEE- Wiley, Piscataway, 
NJ, 1994). 

43 Z. Zhu, S. L. Dvorak, D. L. Heckmann, and J. L. Prince, Radio Sci. 40, 1 (2005). 

44 M. M. Agrest and M. S. Maksimov, Theory of Incomplete Cylindrical Functions and Their 
Applications (Springer, New York, 1971). 

45 S. L. Dvorak, IEEE Antennas Propagat. Mag. 36, 26 (1994). 

46 G. A. Baker Jr. and P. Graves-Morris, Pade Approximants (Cambridge University Press, Cam- 
bridge, 1996). 

47 A. Alu, M. G. Silveirinha, A. Salandrino, and N. Engheta, Phys. Rev. B 75, 155410 (2007). 

48 C. Rizza, A. Di Falco, and A. Ciattoni, Appl. Phys. Lett. 99, 221107 (2011). 

49 H. N. S. Krishnamoorthy, Z. Jacob, E. Narimanov, I. Kretzschmar, and V. M. Menon, Science 
336, 205 (2012). 

50 S. L. Dvorak, IEEE Trans. Microwave Theory Tech. 42, 2164 (1994). 

51 S. L. Dvorak and D. G. Dudley, IEEE Trans. Electromagn. Compat. 37, 192 (1995). 

52 M. Wu, R. G. Olsen, and S. W. Plate, J. Electromagn. Waves Appl. 4, 479 (1990). 

53 E. A. Coddington and N. Levinson, Theory of Ordinary Differential Equations (McGraw-Hill, 
New York, 1955). 

54 Note that the Hq 1 ^ terms in (jBlObj) cancel out when substituted into (1B8|) . 



20 




FIG. 1. (Color online) Problem schematic: a magnetic-current line source is placed at a distance 
x s from a uniaxially-anisotropic metamaterial slab of thickness d and permittivity tensor e given in 
(PQ), immersed in vacuum. Also shown are the global (x,y) and rotated (£, v) Cartesian coordinate 
reference systems, as well as the wavevector k of a propagating plane wave. 



K (a) 
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FIG. 2. (Color online) Typical elliptic (a) and hyperbolic (b) Equi-frequency contours (EFCs) 
pertaining to the dispersion relation in fl5J), for > e v > and < 0,e v > 0, respectively, in the 
global (k x ,k y ) and rotated (k^,k v ) spectral reference systems. 
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FIG. 3. (Color online) Normalized Green's function intensity maps (within and beyond the image 
plane i = d), for a source located x s = Ao/50 away from a slab of thickness d = 0.5Ao, s v = 0, 
eg = —2, and the optical-axis angle a varying from (a) to 75° (f) with step of 15°, computed via 
the reference solution [cf. [5]. Here and henceforth, intensities are normalized with respect to the 
peak- intensity at the image plane in the absence of the slab, i.e., G (d,0,x s ) [cf. Q]. 
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FIG. 4. (Color online) As in Fig. [3j but lateral cuts at the image plane x = d (black-solid curves) 
compared with the "exact" and small-argument-approximated CILHI-based analytical solutions 
(red-dashed and blue-dotted curves, respectively). The magenta-dashed-dotted vertical lines indi- 
cate the fiducial positions y a = dtana of the laterally-displaced images. 




FIG. 5. (Color online) As in Fig. [H but orthogonal cuts at the fiducial position of the image 
y a = dtana. 
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FIG. 6. (Color online) Parameters as in Fig. HI FWHM (empty markers, left axis) and peak- 
intensity (full markers, right axis) at the image plane x = d, estimated via the reference solution 
[cf. ©, circular markers] and the small-argument-approximated CILHI-based analytical solution 
(square markers), as a function of the optical-axis angle a. 
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FIG. 8. (Color online) FWHM (empty markers, left axis) and peak-intensity (full markers, right 
axis) at the image plane x = d, estimated via the reference solution [cf. flSJ)], for x s = Ao/50, 
d=0.5\o, a = 0, Re(e^) = —2, Re(e„) = 10 -3 , as a function of lm(s^ >v ) (in log-scale). Also shown 
(as horizontal dashed lines) are the corresponding estimates from the small-argument-approximated 
CILHI-based analytical solution. 
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FIG. 9. (Color online) As in Fig. [8l but for Re(eg) = —5. Note also the log-scale on the right axis. 
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